Assignment 2¶

by Thomas Campbell¶

Problem 1a¶

How many types of measurements are included in the observation file for each of the GPS, GLONASS, GALIEO and BeiDou systems?

Answer¶

From file:

G 18 C1C L1C D1C S1C C1W S1W C2W L2W D2W S2W C2L L2L D2L S2L C5Q L5Q D5Q S5Q

E 20 C1C L1C D1C S1C C6C L6C D6C S6C C5Q L5Q D5Q S5Q C7Q L7Q D7Q S7Q C8Q L8Q D8Q S8Q

R 20 C1C L1C D1C S1C C1P L1P D1P S1P C2P L2P D2P S2P C2C L2C D2C S2C C3Q L3Q D3Q S3Q

C 12 C2I L2I D2I S2I C7I L7I D7I S7I C6I L6I D6I S6I

G => GPS: 18 measurements.
E => Galileo: 20 measurements.
R => GLONASS: 20 measurements.
C => BeiDou: 12 measurements.

1b.)¶

What are they?

Answer¶

GPS
       L1: {Pseudorange (PR) C/A, Phase C/A, C/A Doppler, C/A Signal Strength, Z-tracking PR, Z-tracking Signal Strength},
       L2: {Z-tracking PR, Z-tracking Phase, Z-tracking Doppler, Z-tracking Signal Strength, L-code PR, L-code Phase, L-code Doppler, L-code Signal Strength},
       L5: {Q PR, Q Phase, Q Doppler, Q Signal Strength}

Galileo
       E1: {C PR, C Phase, C Doppler, C Signal Strength},
       E6: {C PR, C Phase, C Doppler, C Signal Strength},
       E5: {Q PR, Q Phase, Q Doppler, Q Signal Strength},
       E7: {Q PR, Q Phase, Q Doppler, Q Signal Strength},
       E8: {Q PR, Q Phase, Q Doppler, Q Signal Strength}

GLONASS
       G1: {C/A PR, C/A Phase, C/A Doppler, C/A Signal Strength, P PR, P Phase, P Doppler, P Signal Strength},
       G2: {P PR, P Phase, P Doppler, P Signal Strength, C/A PR, C/A Phase, C/A Doppler, C/A Signal Strength},
       G3: {Q PR, Q Phase, Q Doppler, Q Signal Strength}

BeiDou
       B1-2: {I PR, I Phase, I Doppler, I Signal Strength},
       B2b: {I PR, I Phase, I Doppler, I Signal Strength},
       B3: {I PR, I Phase, I Doppler, I Signal Strength}

1c.)¶

Select one satellite from each of the systems (e.g., GPS01, Galileo02, GLONASS15), and plot the corresponding pseudorange measurements (you can choose any code type) as a function of time on a same figure. Make sure the all elements in the figure are clearly labeled

In [ ]:
using DataFrames, Dates, Plots, Statistics
default(seriestype = :scatter, markerstrokewidth=0, xrotation = 45, markersize = .5, right_margin=15Plots.mm)
c = 299792458
F_0 = 10.23
F_1 = 154
F_2 = 120
λ_1 = c / (F_0 * F_1 * 1e6)
λ_2 = c / (F_0 * F_2 * 1e6)


function try_parse(str)
    try
        parse(Float64, str)
    catch e
        return missing
    end
end

function parse_rinex(in_path)

    # read in sp3 file
    in_stream = open(in_path, "r")
    in_lines = readlines(in_stream)

    # the first 22 lines are not needed for these purposes.
    data = in_lines[55:length(in_lines)-1]

    # initialize vectors for each kind of data. sp_d is an intermediate vector, used only to ultimately push to dates.
    dates = []
    entries = []
    sp_d = []
    obs = []

    # iterate over the sp3 dataset. at each row, split, then push to either dates, positions, or velocities vector. 
    # at each occurrence of a position data line, sp_d will also be pushed to dates. 
    # this way, if the number of PRN's changes over the course of the day, the number of repeated dates will change accordingly.
    for line in eachindex(data)
        ln = data[line]

        ln = rpad(ln, 321, " ")

        sp = split(ln)
        if sp[1] == ">"
            sp_d = copy(sp)

        else
            push!(entries, try_parse.([ln[4:17], ln[20:33], ln[36:49], ln[52:65], ln[68:81], ln[84:97], ln[100:113],
                ln[116:129], ln[132:145], ln[148:161], ln[164:177], ln[180:193], ln[196:209],
                ln[212:225], ln[228:241], ln[244:257], ln[260:273], ln[276:289],
                ln[292:305], ln[308:321]]))

            push!(dates, parse.(Float64, sp_d[2:end-2]))
            push!(obs, ln[1:3])
        end
    end

    # create a matrix of parsed sp3 data
    return vcat.(dates, obs, entries)

end

function make_rindf(rinex, gnss)

    col_names = Dict("gps" => ["C1C", "L1C", "D1C", "S1C", "C1W", "S1W", "C2W", "L2W", "D2W", "S2W", "C2L", "L2L", "D2L", "S2L", "C5Q", "L5Q", "D5Q", "S5Q"],
        "galileo" => ["C1C", "L1C", "D1C", "S1C", "C6C", "L6C", "D6C", "S6C", "C5Q", "L5Q", "D5Q", "S5Q", "C7Q", "L7Q", "D7Q", "S7Q", "C8Q", "L8Q", "D8Q", "S8Q"],
        "sbas" => ["C1C", "L1C", "D1C", "S1C", "C5I", "L5I", "D5I", "S5I"],
        "glonass" => ["C1C", "L1C", "D1C", "S1C", "C1P", "L1P", "D1P", "S1P", "C2P", "L2P", "D2P", "S2P", "C2C", "L2C", "D2C", "S2C", "C3Q", "L3Q", "D3Q", "S3Q"],
        "bds" => ["C2I", "L2I", "D2I", "S2I", "C7I", "L7I", "D7I", "S7I", "C6I", "L6I", "D6I", "S6I"],
        "qzss" => ["C1C", "L1C", "D1C", "S1C", "C2L", "L2L", "D2L", "S2L", "C5Q", "L5Q", "D5Q", "S5Q"],
        "irnss" => ["C5A", "L5A", "D5A", "S5A"])

    if gnss == "gps"
        gnss_flag = 'G'
    elseif gnss == "galileo"
        gnss_flag = 'E'
    elseif gnss == "sbas"
        gnss_flag = 'S'
    elseif gnss == "glonass"
        gnss_flag = 'R'
    elseif gnss == "bds"
        gnss_flag = 'C'
    elseif gnss == "qzss"
        gnss_flag = 'J'
    elseif gnss == "irnss"
        gnss_flag = 'I'
    end

    gnss_data = []

    for line in eachindex(rinex)
        if rinex[line][7][1] == gnss_flag
            push!(gnss_data, rinex[line])
        end
    end

    full_names = ["year", "month", "day", "hour", "minute", "seconds", "sat"]
    date_names = copy(full_names[1:end-1])
    append!(full_names, col_names[gnss])

    df = DataFrame(mapreduce(permutedims, vcat, gnss_data), :auto)
    df = df[:, collect(1:length(full_names))]

    rename!(df, Symbol.(full_names))

    df.dt = DateTime.(Dates.Year.(df.year),
        Dates.Month.(df.month),
        Dates.Day.(df.day),
        Dates.Hour.(df.hour),
        Dates.Minute.(df.minute))

    select!(df, Not(date_names))
    return df
end

make_rindf (generic function with 1 method)
In [ ]:
in_path = "C:/Users/tom_j/Downloads/ABMF00GLP_R_20200200000_01D_30S_MO.rnx"

# Parse Rinex file
parsed = parse_rinex(in_path);

# Create Datasets
mygps = make_rindf(parsed, "gps");
mygalileo = make_rindf(parsed, "galileo");
myglonass = make_rindf(parsed, "glonass");
mybeidou = make_rindf(parsed, "bds");

# Pick a satellite
g30 = dropmissing(mygps[mygps.sat .== "G30", :]);
e11 = mygalileo[mygalileo.sat .== "E11", :];
r12 = myglonass[myglonass.sat .== "R12", :];
c35 = mybeidou[mybeidou.sat .== "C35", :];
In [ ]:
plot(g30.dt, g30.C1C, label= "GPS L1 C/A", ylabel = "Range (m)", title = "GNSS Pseudoranges")
plot!(e11.dt, e11.C1C, label="Galileo E1 C/A")
plot!(r12.dt, r12.C1C,label="GLONASS G1 C/A")
plot!(c35.dt, c35.C2I, label="BeiDou B1-2 I")

Problem 2a.)¶

Plot phase L1 and code P1 on a same plot for each of the satellite in problem 1.

In [ ]:
plot(g30.dt, g30.L1C * λ_1, label="GPS L1 Phase", ylabel="Range (m)", title="GPS P Code PR vs L1 Phase")
plot!(g30.dt, g30.C1W, label="GPS P1 Pseudorange")
plot!(twinx(),g30.dt, g30.L1C * λ_1 - g30.C1W , label= "Φ - PR", legend=:topleft, markercolor = :black)
In [ ]:
plot(e11.dt, e11.L1C * λ_1,label="Galileo E1 C Phase", ylabel="Range (m)", title="Galileo E5 Code PR vs E1 Phase", legend=:topleft)
plot!(e11.dt, e11.C5Q,label="Galileo E5 Q Pseudorange")
plot!(twinx(), e11.dt, e11.L1C * λ_1 - e11.C5Q, label="Phase - PR",legend=:bottomleft, markercolor = :black)
In [ ]:
plot(r12.dt, r12.L1C * c / (1602*1e6),label="GLONASS G1 C Phase", ylabel="Range (m)", title="GLONASS G2a Code PR vs G1 Phase")
plot!(r12.dt, r12.C2P, label="GLONASS G2a P Pseudorange")
plot!(twinx(), r12.dt, r12.L1C * c / (1602*1e6) - r12.C2P, label="Phase - PR", legend=:bottomleft, markercolor=:black)
In [ ]:
plot(c35.dt, c35.L2I * c / (1561.098 * 1e6), label="BeiDou B1-2 I Phase", ylabel="Range (m)", title="BeiDou C6I Code PR vs B1-2 Phase")
plot!(c35.dt, c35.C6I, label="BeiDou B3 I Pseudorange")
plot!(twinx(), c35.dt, c35.L2I * c / (1561.098 * 1e6) - c35.C6I, label="Phase - PR", legend=:bottomleft, markercolor=:black)

2b.)¶

Would you conclude one is more precise than the other? And why?

Answer¶

It seems hard to tell which measurement would be more precise from the plots above.

Need to look at the chip rates and precision for carrier phases. Maybe they are close given the different codes?

GNSS Constellation or Measurement?

Problem 3a.)¶

Why can ionosphere-free combinations remove the first-order ionospheric effect?

3b.)¶

Form code (Pc ) and carrier phase (L c ) ionosphere-free combinations for the three satellites.

In [ ]:
function ion_free(f_1, f_2, m_1, m_2)
    γ = f_1/f_2
    α = γ*γ/(1 - γ*γ)
    return m_1 + 2*α*(m_1 - m_2)
end

#plot(g30.dt, ion_free(F_1, F_2, g30.C1C, g30.C2W), markerstrokewidth=0, xrotation=45, label="Code IF", ylabel="Range (m)", legend=:topleft)
plot(g30.dt, ion_free(F_1, F_2, λ_1 * g30.L1C, λ_2 * g30.L2W), label="Phase IF")
plot!(g30.dt, λ_1 *g30.L1C,label="GPS P1 Pseudorange")

Problem 4a.)¶

The code multipath can be check by plotting the difference of code and carrier ionosphere-free combinations (Pc - L c ). Form the combination for the three satellites.

In [ ]:
plot(g30.dt, ion_free(F_1, F_2,g30.C1C, g30.C1W) - ion_free(F_1, F_2, λ_1 * g30.L1C, λ_2 * g30.L2W), label="Phase IF")

4b.)¶

What can you tell from the results?

Problem 5a.)¶

What is geometry-free combination, and what is the physical meaning of this combination?

5b.)¶

What are the main utilities of the combination?

5c.)¶

Form the geometry-free combination, e.g., L1-L2, P1-P2 and P1-L1, for the three satellites.

In [ ]:
plot(g30.dt, g30.C1C - g30.C1W)
plot!(g30.dt, λ_1 * g30.L1C - λ_2 * g30.L2W )
plot!(g30.dt, g30.C1W - λ_1 * g30.L1C)

#Why are these all shifted around?

5d.)¶

Among three combinations (e.g., L1-L2, P1-P2 and P1-L1), which one shows a greater level of noise, and why?

5e.)¶

How would you explain the variations in these combinations with time?

5f.)¶

Why do the combinations of L1-L2 and P1-P2 have opposite sign?

5g.)¶

Plot the combination (L1-L2) and (P2-P1) on the same plot.

In [ ]:
plot(g30.dt, g30.C1C - g30.C1W)
plot!(g30.dt,λ_2 * g30.L2W - λ_1 * g30.L1C )

5h.)¶

Does dispersion increase at the beginning and end of arcs?

Yes

5i.)¶

How would you explain the increase if you do see the increase?

Arcs usually begin and end at low elevations. At low elevations, signals are more likely to experience obstructions, increased tropospheric and ionospheric delays, as well as more multipath.